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Mixing in a stratified shear flow: 
energetics and sampling 

By G. N. Ivey 1 , J. R. Koseff 2 , D. A. Briggs 2 AND J. H. Ferziger 2 

Direct numerical simulations of the time evolution of homogeneous stably strat- 
ified shear flows have been performed by Holt (1990), Holt et al. (1992), and 
Itsweire et al. (1992) for Richardson numbers from 0 to 1 and for Prandtl numbers 
between 0.1 and 2. The results indicate that mixing efficiency Rj varies with turbu- 
lent Froude number in a manner consistent with laboratory experiments performed 
with Prandtl numbers of 0.7 and 700. However, unlike the laboratory results, for 
a particular Froude number, the simulations do not show a clear dependence on 
the magnitude of Rj on Pr. The observed maximum value of Rf is 0.25. When 
averaged over vertical length scales of an order of magnitude greater than either 
the overturning or Ozmidov scales of the flow, the simulations indicate that the 
dissipation rate e is only weakly lognormally distributed with an intermittency of 
about 0.01 whereas estimated values in the ocean (Baker and Gibson (1987)) are 3 
to 7. 

1. Introduction 

The specification of the buoyancy flux in a fluid with a stable density gradient 
and subjected to forcing by either mean shear or a mechanical means of generating 
turbulent kinetic energy is crucial to the understanding of mixing processes in a 
wide variety of geophysical applications. Quantifying the rate of vertical mixing in 
the ocean, for example, is central not only to parameterizing vertical heat and salt 
fluxes, but also those of passive tracers such as pollutants and chemicals. Thus, 
predicting the mixing efficiency or flux Richardson number Rj is of central im- 
portance in turbulent mixing in density stratified flows. Ivey and Imberger (1991) 
(hereafter II) introduced a generalized definition of Rf based on the full turbulent 
kinetic energy equation and demonstrated that available laboratory measurements 
had peak values of Rf = 0.20 for fluids with Prandtl number of 700 and 0.15 for 
fluids with Prandtl number of 0.7. In the laboratory work, turbulence was gener- 
ated with grids and, in the case of Rohr et al. (1988), by both grids and a mean 
shear. This work has been extended by the direct numerical simulation work of 
Holt (1990), Holt et al. (1992), and Itsweire et al. (1992), who numerically studied 
the temporal evolution of homogeneous turbulence subjected to a constant mean 
velocity and density gradients using the pseudo-spectral method of Rogallo (1981). 
The Boussinesq form of the Navier Stokes equations were solved for the three di- 
mensional velocity and density fields on a 128 x 128 x 128 grid. The velocity fields 
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were initially isotropic with Reynolds numbers based on the integral scale and the 
turbulent kinetic energy ranging from 52 to 104, Richardson numbers ranging from 
0.058 to 1.0, and Prandtl numbers of either 1 or 2. In the oceanic case, the density 
stratification is predominantly due to vertical temperature gradients. Therefore, in 
order to investigate the differences in the behavior of Rj with Pr observed by II, 
we have extended the earlier work of Holt et al (1992) by performing additional 
simulations at values of Ri from 0.075 to 1.0 with Prandtl numbers as low as 0.1; 
these results are described in Section 3 below. 

Oceanic turbulence is intermittent in character, and there exists considerable con- 
troversy about how to correctly sample the flow to determine the volume averaged 
dissipation (eg. Baker and Gibson (1987)). We examine this question by averaging 
over sub- sections of our 128 x 128 x 128 domain and calculating the intermit tency, 
defined as the variance of the logarithm of a quantity: the results are summarized 
in Section 4. Finally, Section 5 summarizes our conclusions. 


2. Dynamics of mixing in stratified flows 

The full turbulent kinetic energy equation can be written as: 


dt { 2 ’ dxj 2 



1 , dUi 
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where q 2 = u 2 + u\ + u\ is twice the turbulent kinetic energy, p Q is the background 
density, e = is the rate of dissipation of turbulent kinetic energy, and 

velocity component is in the direction of the gravity vector. For homogeneous 

turbulence at steady state, (1) can be written as P = B+e where P is the production 
of turbulent kinetic energy (the last term on the left hand side of (1)) Then, II 
defined 

(2) 

and also introduced the concept of a turbulent Froude number associated with the 
energy bearing eddies defined by 


Fry = 


sr 
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where L c is the scale of the most energetic overturns, L 0 = ( e/N 3 ) represents 

the Ozmidov scale, and N is the buoyancy frequency characterizing the stable back- 
ground stratification. An alternate form of (2) is 


R f 


1 

1 + Fr 2 T jR pw 


(4) 


where R pw is the correlation coefficient between the density and vertical velocity 
fluctuations. II found the available laboratory measurements were in accord with 
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Froude number, Frp 

Figure 1 . Rf as a function of Fr T = ( L a /L e ) 2/ 3 for simulations with Pr = 0.1 
(o), Pr = 0.5 (o) and Pr = 2.0 (•). Experimental data is provided by Stillinger et 
al., Itsweire et al ., and Rohr et al. for Pr = 700 (x) and Lienhard and Van Atta 
for Pr = 0.7 (+). 


the simple predictions in (4). Additionally, II introduced two further dimensionless 
numbers, the turbulent Reynolds number Rex and small scale Froude number Fr 7 


defined by 



i? er = UjLc =(f c ) 4 / 3 

v Lk 

(5) 

and 


(6) 


where e/vN 2 is a measure of the range of overturning scales when buoyancy strongly 
affects the flow and L* = {u 3 represents the Kolmogorov scale. The dimen- 

sionless numbers (3), (5), and (6) are related by 


Fry = ( 


1 

ReT 


) 1/2 Fr y . 


(7) 


One final length scale is used in discussing the numerical results below. The Ellison 


scale L t is defined as 


L e 





( 8 ) 


where ~p{z) is mean density. As Itsweire et al. (1992) show, except at very high 
Richardson numbers (based on mean shear and density gradient), L c and L e are 
the same scale to within a constant of 0(1). 
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FIGURE 2. as a function of Frx for Pr = 2. 

3. Mixing efficiency 

All results in this section were computed by ensemble- averaging data over the 
entire computational domain at each time step. In order to eliminate transient 
behavior associated with the initial conditions, only the data for St > 2 are con- 
sidered, where the shear time, St, is a dimensionless time with S = dUjdz . Rf 
was computed according to the definition in (2), as were the dimensionless num- 
bers Frx, Rer , and JFY 7 from their respective definitions in equations (3), (5), and 
(6). Figure 1 summarizes the results for the simulations with Pr = 0.1, 0.5, and 
2. Mixing efficiency, i?/, is plotted against Frr to facilitate direct comparison with 
the laboratory observations (presented by II) which are also shown in Figure 1. 
The numerical results show the same general distribution as the laboratory mea- 
surements. As Frr becomes very large, the overturning scale is small compared 
to the Ozmidov scale, the mixing becomes inefficient because there is much more 
turbulent kinetic energy than necessary to mix the scalar gradient, and Rf rapidly 
decreases. Similarly, for very small values of Fr t when buoyancy effects dominate 
and reduce the mixing, Rf again sharply decreases. At an intermediate value of 
Frr in the range of 1 to 1.5, the mixing is optimal with peak values of Rf = 0.25. 

The behavior of R pw as a function of Frr given in Figure 2 for Pr = 2 (as 
an example) shows a similar trend. For Frx > 1, the correlation coefficient R pw 
tends to a constant value of about 0.6; hence, from Equation (4), Rf will rapidly 
decrease with increasing Frr . Conversely, for Frr < 1, R pw decreases rapidly as 
the density and velocity fluctuations become uncorrelated, and Rf tends to zero as 
a consequence of the denominator in (4) being much greater than unity. 

While the overall behavior in the numerical and laboratory results is the same 
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Figure 3. R / as a function of c/uN 2 for simulations with Pr = 0.1 (o), 0.5 (•), 
2.0 (x). 

in Figure 1, there are differences in detail. Unlike the laboratory data, for a given 
Frr, the numerically derived values of Rf are independent of Prandtl number and 
do not show the experimental tendency of Rf to increase with Pr. Nevertheless, 
recalling that the laboratory data are derived from time-averaged statistics in a 
steady mean flow whereas the numerical value are ensemble- averaged values in an 
evolving flow, the differences are remarkably small. The implication is that the 
peak mixing efficiency is 0.25 for Fr r 1 to 1.5, irrespective of Prandtl number. 

In field measurements of oceanic turbulence, the overturning scales L t or L c are 
not usually measured while e invariably is. Substituting (7) into (4) yields the more 
practical form 

Rf = 1 + Wf^N*) 

where (3 — R~^Rex- 

In Figure 3, we re-plot the results from the simulations against e/uN 2 (= Fr 2 ). 
The minimum dissipation needed to sustain a vertical buoyancy flux, and hence 
positive i?/, is clearly a strong function of Pr. For Pr = 0.1, e/uN 2 may be as 
little as 0.4 and still sustain a positive buoyancy flux , whereas for Pr = 2, the 
minimum e/uN 2 for a positive buoyancy flux is about 20. The other interesting 
feature of Figure 3 is that Rf becomes independent of Pr for large ejvN 2 . For 
e/uN 2 > 10, a best fit is 0 = (e/vN 2 y°' 6 , hence (9) simplifies to 


1 

f ~l + (e/uN 2 ) 0A 


( 10 ) 
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dissipation, ecm 2 s 3 

FIGURE 4. Lognormal plot of dissipation estimates with no averaging (Ri = .075, 
Pr = 2, St — 6). 

and using (2) 

B = e°- 6 (uN 2 )° 4 (11) 

which provides an approximate but simple means of computing buoyancy flux from 
the dissipation for relatively energetic flows. 

4. Sampling turbulence in a stratified fluid 

Turbulence measurements are made in the ocean with either vertically falling mi- 
crostructure instruments or, less commonly, horizontally towed instruments. The 
buoyancy flux B is not directly measured but, as indicated above, dissipation esti- 
mates are made and then B is computed by estimating Rf as outlined in Section 3 
(see also Itsweire et ah 1992). For falling probes, dissipation estimates are typically 
made by measuring two turbulent velocity components and computing total dissi- 
pation using models (see Itsweire et ah 1992). This procedure produces estimates 
of dissipation averaged over about 2 meters in the vertical, and these estimates are 
further averaged to characterize the dissipation on much larger scales such as the 
oceanic thermocline (for example see review of Gregg 1987). Gibson and Baker 
(1987) and Gibson (1991) have argued that dissipation in oceanic turbulence is log- 
normally distributed with an intermittency, er 2 , in the range of 3 to 7. Furthermore, 
they claim that, due to the large scales and the limited sampling, the dissipation is 
greatly undersampled. Given the large intermittency, Baker and Gibson maintain 
that to estimate of e to within ±10% one would need to average the dissipation 
calculated from thousands of independent sampling profiles! Gurvich and Yaglom 
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(1967) (see also Yamazaki and Lueck (1990)) developed three criteria which must 
be satisfied in order for dissipation to be a lognormally distributed quantity in the 
ocean. 

(i) The turbulence must be homogeneous, 

(ii) the averaging scale, r, must be small compared to the domain scale, L, and 

(iii) the averaging scale must be large compared to the Kolmogorov scale Lk- 

Yamazaki and Lueck (1990) suggest that (iii) can be satisfied with r as low as 
3 Z*, but in all oceanic datasets, there is uncertainty about meeting criterion (i) 
given the patchy nature of the turbulence away from regions of highly energetic 
forcing such as the near-surface region. We investigated the sampling question by 
analyzing the results from one typical simulation. In particular, we examined the 
numerical results corresponding to Ri = 0.075, Pr = 2 at St = 6. For this data 
set, L = 25 cm, the grid scale, r = Lj 128, is 0.195 cm, the ensemble- averaged 
dissipation is 0.213 cm 2 s~ 3 , Lk = 0.0465 cm, L c = 2.68 cm, and L e = 1.37 cm. For 
the full data set of 128 3 points without averaging, all three criteria of Yamazaki and 
Lueck (1990) are met since r = 4.2 Lk- The corresponding distribution of dissipation 
should, therefore, be lognormal, and in Figure 4 the strong lognormality of the data 
with a 2 ~ .75 is evident. However, for comparison with typical oceanic sampling, 
of greater interest are the consequences of averaging estimates in the vertical. We 
chose to examine the effect of averaging over the full depth of the computational 
domain (25 cm), which is equal to 18.2 L e or 9.3 L 0 . Even though we are averaging 
over the depth, we still have a statistically significant sample of 128 2 points. Note 
that we do not satisfy criterion (ii) because r = L. 

Figure 5 shows the distribution of vertically averaged e for four different Richard- 
son numbers. The effect of averaging is quite dramatic and the value of a 2 is reduced 
to 0.01 and becomes essentially independent of Ri. Thus, when averaged over ver- 
tical scales significantly greater than L e or L c , dissipation is only weakly lognormal 
with cr 2 less than the value of 3 to 7 suggested for the ocean. Using techniques 
developed in Baker and Gibson, this translates to about 4 to 8 required profiles to 
obtain an estimate of the mean value of e within ±10%. 

In order to more fully explore the second criterion of Gurvich and Yaglom, we 
averaged the data in the vertical but over length scales r which were smaller than 
L, e.g. 16 points or about 2.5 L c . The distribution of e averaged over this scale is 
plotted in Figure 6 which shows that the data is indeed lognormal but with a 2 ~ .1, 
much smaller than the 3 - 7 obtained by Baker and Gibson (1987). Again, the 
implication is that, if the data are averaged over length scales comparable to the 
overturning scale, the intermittency will not be as high as that described by Baker 
and Gibson (1987). This implies that much fewer sampling profiles are necessary 
to obtain a reasonable estimate of e. Finally, we should note that since oceanic 
samples are likely to be drawn from a domain with significant regions of minimal 
dissipation the homogeneity criteria may be the most difficult restriction to satisfy. 
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FIGURE 5. Lognormal probability plot of e for Ri = 0.075 (o), 0.21 (o), 0.37 (x) 
and 1.0 (+) with Pr = 2.0, St = 6. 
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FIGURE 6. Lognormal probability plot of ln(e) with the averaging scale r ~ 2.5L e , 
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5. Conclusions 

Direct numerical simulations of homogeneous turbulence in stably stratified shear 
flows confirm the variation of flux Richardson number Rf with turbulent Froude 
number Frr and e/vN 2 observed in laboratory experiments. With Rf defined 
as the buoyancy flux divided by the production of turbulent kinetic energy, there 
appears to be no systematic dependence of Rf on Pr in the range of Pr from 0.1 
to 2. This result is not consistent with the laboratory observations; however, the 
differences in Rf between the simulations and the experiments are small, and the 
data from all sources indicate that Rf has a peak of 0.25, independent of Pr. Sub- 
sampling of the computational domain of 128 3 points was investigated to examine 
the distribution of the dissipation. The results indicate that when dissipation is 
estimated by averaging over vertical scales of an order of magnitude greater than 
either the Ellison or Ozmidov scales, the distribution is very weakly lognormal with 
an intermittency, a 2 ~ 0.01. This value is considerably smaller than some estimates 
in the oceanic literature and suggests sampling restrictions may not be as severe 
as previously suggested provided the sampling and averaging are performed over 
domains where the turbulence is homogeneous. 

Acknowledgements 

The authors are very grateful to the CTR for making this work possible. JRK, 
DAB and JHF also wish to acknowledge the Office of Naval Research for their 
support of this work through grant number N00014-92-J-1611. 

REFERENCES 

BAKER, M. & Gibson, C. G. 1987 Sampling turbulence in the stratified ocean: 
statistical consequences of strong intermittency. J. Phys. Oceanogr. 17, 417- 
440. 

GlBSON, C. H. 1991 Turbulence and Stochastic Processes: Kolmogorov’s ideas 50 
years on, Hunt, J. C., Phillips, O. M. k Williams, D (eds). The Royal Society, 
London. 

Gregg, M. C. 1987 Diapycnal mixing in the thermocline: A review. J. Geophys. 
Res. 92, 5249-5286. 

GURVICH, A. S. k YagLOM, A. M. 1967 Breakdown of eddies and probability 
distributions for small scale turbulence. Phys. Fluids. 10, 59-65. 

Holt, S. E., Koseff, J. R. &: FERZIGER, J. H. 1992 The evolution of turbulence 
in the presence of mean shear and stable stratification. J. Fluid Mech. 237 , 499- 
539. 

LlENHARD, J. H. k VAN Atta, C. W. 1990 The decay of turbulence in thermally 
stratified flow. J. Fluid Mech. 210, 57-112. 

Itsweire, E. C., Helland, K. N. k Van Atta, C. W. 1986 The evolution 
of grid-generated turbulence in a stably stratified flow. J. Fluid Mech. 162, 
299-338. 



344 


G. N. Ivey, J. R. Koseff, D. A. Briggs & J. H. Ferziger 


Itsweire, E. C., Koseff, J. R., Briggs, D. A. Sc Ferziger, J. H. 1992 

Turbulence in stratified shear flows: Implications for interpreting shear-induced 
mixing in the ocean. J. Phys. Oceanogr, ( accepted for publication). 

IVEY, G. N. Sc Imberger, J. 1991 On the nature of turbulence in a stratified 
fluid. Part I: Energetics of mixing. J. Phys. Oceanogr. 21, 650-658. 

ROGALLO, R. S. 1981 Numerical experiments in homogeneous turbulence. NASA 
Tech. Memo 81315. 

Rohr, J. J., Itsweire, E. C., Helland, K. N. Sc Van Atta, C. W. 1988 An 
investigation of the growth of turbulence in a uniform-mean-shear flow. J. Fluid 
Mech. 187, 1-33. 

Stillinger, D. C., Helland, K. N. Sc Van Atta, C. W. 1983 Experiments on 
the transition of homogeneous turbulence to internal waves in a stratified fluid. 
J. Fluid Mech. 131, 91-122. 

YAMAZAKI, H. Sc Lueck, P. 1990 Why oceanic dissipation rates are not lognor- 
mal. J. Phys. Oceanogr. 20(12), 1907-1918. 


